
#include "SinglePyrolysisSolidADMaterial.h"

registerMooseObject("TrilobitaApp", SinglePyrolysisSolidADMaterial);

InputParameters
SinglePyrolysisSolidADMaterial::validParams()
{
    auto params = ADMaterial::validParams();

    params.addRequiredCoupledVar("temperature", "Coupled temperature");
 //   params.addRequiredCoupledVar("rhos", "Coupled solid density");
    params.addRequiredParam<Real>("rhov", "density of original material");
    params.addRequiredParam<Real>("rhoc", "density of fully charred material");
    params.addRequiredParam<Real>("T_begin", "temperature at which the reaction begins");
    params.addRequiredParam<Real>("precofficient", "pre-exponential factor");
    params.addRequiredParam<Real>("ER", "E/R");
    params.addRequiredParam<Real>("reaction_level", "Reaction level of pyrolysis,normally 1,2 or 3");
    params.addRequiredParam<std::vector<Real> >("kv_list", "The vector of kv values of origin material for building the piecewise function");
    params.addRequiredParam<std::vector<Real> >("kvT_list", "The vector of temperature values for building the piecewise function to compute kv ");
    params.addRequiredParam<std::vector<Real> >("kc_list", "The vector of kc values for building the piecewise function to compute kc ");
    params.addRequiredParam<std::vector<Real> >("kcT_list", "The vector of temperature values for building the piecewise function to compute kc ");
    params.addRequiredParam<std::vector<Real> >("cpv_list", "The vector of cpv values for building the piecewise function to compute thermal conductivity of origin material ");
    params.addRequiredParam<std::vector<Real> >("cpvT_list", "The vector of cpv values for building the piecewise function to compute thermal conductivity of origin material ");
    params.addRequiredParam<std::vector<Real> >("cpc_list", "The vector of cpc values for building the piecewise function to compute cpc ");
    params.addRequiredParam<std::vector<Real> >("cpcT_list", "The vector of temperature values for building the piecewise function to compute cpc");
    params.addRequiredParam<Real>("pyrolysis_heat", "pyrolysis heat");
    params.addRequiredParam<Real>("porosity_v", "porosity of original material");
    params.addRequiredParam<Real>("porosity_c", "porosity of fully charred material");
    params.addRequiredParam<Real>("permeability_v", "permeability of original material");
    params.addRequiredParam<Real>("permeability_c", "permeability of fully charred material");
    params.addRequiredParam<Real>("pyrolysis_heat", "pyrolysis heat");
    params.addRequiredParam<Real>("surface_emissivity", "emissivity");

    return params;
}

SinglePyrolysisSolidADMaterial::SinglePyrolysisSolidADMaterial(const InputParameters& parameters)
    : ADMaterial(parameters),
    _rhos(declareADProperty<Real>("solid_density")),
    _pyrolysis_source(declareADProperty<Real>("pyrolysis_source")),
    _drhos_dt(declareADProperty<Real>("drhos_dt")),
    _rhov(declareADProperty<Real>("rhov")),
    _rhoc(declareADProperty<Real>("rhoc")),
    _char_percent(declareADProperty<Real>("char_percent")),
    _thermal_conductivity_v(declareADProperty<Real>("solid_thermal_conductivity_v")),
    _thermal_conductivity_c(declareADProperty<Real>("solid_thermal_conductivity_c")),
    _thermal_conductivity_s(declareADProperty<Real>("solid_thermal_conductivity")),
    _specific_heat_v(declareADProperty<Real>("solid_specific_heat_v")),
    _specific_heat_c(declareADProperty<Real>("solid_specific_heat_c")),
    _specific_heat_s(declareADProperty<Real>("solid_specific_heat")),
    _porosity_v(declareADProperty<Real>("porosity_v")),
    _porosity_c(declareADProperty<Real>("porosity_c")),
    _porosity(declareADProperty<Real>("porosity")),
    _permeability_v(declareADProperty<Real>("permeability_v")),
    _permeability_c(declareADProperty<Real>("permeability_c")),
    _permeability(declareADProperty<Real>("permeability")),
    _pyrolysis_heat(declareADProperty<Real>("pyrolysis_heat")),
    _surface_emissivity(declareProperty<Real>("surface_emissivity")),
    _rhos_old(getMaterialPropertyOld<Real>("solid_density")),
    _temperature(coupledValue("temperature"))
 //   _rhos_value(coupledValue("rhos"))
{
    _rhov_value = getParam<Real>("rhov");
    _rhoc_value = getParam<Real>("rhoc");
    _T_begin_value = getParam<Real>("T_begin");
    _precoff_value = getParam<Real>("precofficient");
    _ER_value = getParam<Real>("ER");
    _reaction_level_value = getParam<Real>("reaction_level");
    _kvT_list = getParam<std::vector<Real>>("kvT_list");
    _kv_list = getParam<std::vector<Real>>("kv_list");
    _kcT_list = getParam<std::vector<Real>>("kcT_list");
    _kc_list = getParam<std::vector<Real>>("kc_list");
    _cpvT_list = getParam<std::vector<Real>>("cpvT_list");
    _cpv_list = getParam<std::vector<Real>>("cpv_list");
    _cpcT_list = getParam<std::vector<Real>>("cpcT_list");
    _cpc_list = getParam<std::vector<Real>>("cpc_list");
    _func_kv_T = LinearInterpolation(_kvT_list, _kv_list);
    _func_kc_T = LinearInterpolation(_kcT_list, _kc_list);
    _func_cpv_T = LinearInterpolation(_cpvT_list, _cpv_list);
    _func_cpc_T = LinearInterpolation(_cpcT_list, _cpc_list);
    _pyrolysis_heat_value = getParam<Real>("pyrolysis_heat");
    _porosity_v_value = getParam<Real>("porosity_v");
    _porosity_c_value = getParam<Real>("porosity_c");
    _permeability_v_value = getParam<Real>("permeability_v");
    _permeability_c_value = getParam<Real>("permeability_c");
    _pyrolysis_heat_value = getParam<Real>("pyrolysis_heat");
    _surface_emissivity_value = getParam<Real>("surface_emissivity");
}

void
SinglePyrolysisSolidADMaterial::initQpStatefulProperties()
{
   
    _rhos[_qp] = _rhov_value;
    _rhov[_qp] = _rhov_value;
    _rhoc[_qp] = _rhoc_value;
    _thermal_conductivity_v[_qp] = _func_kv_T.sample(_temperature[_qp]);
    _thermal_conductivity_c[_qp] = _func_kc_T.sample(_temperature[_qp]);
    _specific_heat_v[_qp] = _func_cpv_T.sample(_temperature[_qp]);
    _specific_heat_c[_qp] = _func_cpc_T.sample(_temperature[_qp]);
    _char_percent[_qp] = (_rhov_value - _rhos[_qp]) / (_rhov_value - _rhoc_value + 1e-8);
    _thermal_conductivity_s[_qp] = _thermal_conductivity_v[_qp] * (1 - _char_percent[_qp]) + _thermal_conductivity_c[_qp] * _char_percent[_qp];;
    _specific_heat_s[_qp] = _specific_heat_v[_qp] * (1 - _char_percent[_qp]) + _specific_heat_c[_qp] * _char_percent[_qp];
    _porosity_v[_qp] = _porosity_v_value;
    _porosity_c[_qp] = _porosity_c_value;
    _permeability_v[_qp] = _permeability_v_value;
    _permeability_c[_qp] = _permeability_c_value;
    _porosity[_qp] = _porosity_v[_qp] * (1 - _char_percent[_qp]) + _porosity_c[_qp] * _char_percent[_qp];
    _permeability[_qp] = _permeability_v[_qp] * (1 - _char_percent[_qp]) + _permeability_c[_qp] * _char_percent[_qp];
    _pyrolysis_heat[_qp] = _pyrolysis_heat_value;
    _surface_emissivity[_qp] = _surface_emissivity_value;
}
Real SinglePyrolysisSolidADMaterial::PyrolysisSource(Real rho, Real T, Real T_begin, Real rhoc, Real rhov, Real precoff, Real ER, Real reaction_level)
{
    if (T > T_begin)
    {
        if ((rho - rhoc) / (rhov - rhoc) < 0.01)
        {
            return 0;
        }
        else
        {
            Real tmp = (rho - rhoc) / rhov;
            Real r = precoff * exp(-ER / T) * rhov * pow(tmp, reaction_level);
            return precoff * exp(-ER / T) * rhov * pow(tmp, reaction_level);
        }
    }
    else
    {
        return 0;
    }
}

void
SinglePyrolysisSolidADMaterial::computeQpProperties()
{
    _pyrolysis_source[_qp] = -PyrolysisSource(_rhos_old[_qp], _temperature[_qp], _T_begin_value, _rhoc_value, _rhov_value, _precoff_value, _ER_value, _reaction_level_value);
    _drhos_dt[_qp] = _pyrolysis_source[_qp];
    _rhov[_qp] = _rhov_value;
    _rhoc[_qp] = _rhoc_value;
   if( (_rhos_old[_qp] + _drhos_dt[_qp] * _dt) > _rhoc[_qp])
    {
        _rhos[_qp] = _rhos_old[_qp] + _drhos_dt[_qp] * _dt;
    }
    else
    {
        _rhos[_qp] = _rhoc[_qp];
    }
 //   _rhos[_qp] = _rhos_value[_qp];
    _thermal_conductivity_v[_qp] = _func_kv_T.sample(_temperature[_qp]);
    _thermal_conductivity_c[_qp] = _func_kc_T.sample(_temperature[_qp]);
    _specific_heat_v[_qp] = _func_cpv_T.sample(_temperature[_qp]);
    _specific_heat_c[_qp] = _func_cpc_T.sample(_temperature[_qp]);
    _char_percent[_qp] = (_rhov_value - _rhos[_qp]) / (_rhov_value - _rhoc_value + 1e-8);
    _thermal_conductivity_s[_qp] = _thermal_conductivity_v[_qp] * (1 - _char_percent[_qp]) + _thermal_conductivity_c[_qp] * _char_percent[_qp];;
    _specific_heat_s[_qp] = _specific_heat_v[_qp] * (1 - _char_percent[_qp]) + _specific_heat_c[_qp] * _char_percent[_qp];
    _porosity_v[_qp] = _porosity_v_value;
    _porosity_c[_qp] = _porosity_c_value;
    _permeability_v[_qp] = _permeability_v_value;
    _permeability_c[_qp] = _permeability_c_value;
    _porosity[_qp] = _porosity_v[_qp] * (1 - _char_percent[_qp]) + _porosity_c[_qp] * _char_percent[_qp];
    _permeability[_qp] = _permeability_v[_qp] * (1 - _char_percent[_qp]) + _permeability_c[_qp] * _char_percent[_qp];
    _pyrolysis_heat[_qp] = _pyrolysis_heat_value;
    _surface_emissivity[_qp] = _surface_emissivity_value;

}